#include "Equation.h"
#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

double f(double x)
{
	return exp(sin(x)) * (sin(x) - cos(x) * cos(x));
	//return (-x * x + x + 2) * sin(x) + (4 * x - 2) * cos(x);
}

double u(double x)
{
	return exp(sin(x));
	//return x * (1-x) * sin(x);
}

const Vector GetBoundary(std::function<double(double x)> _u)
{
	Vector v(2);
	v[0] = u(0);
	v[1] = u(1);
	return v;
}

double MaxNorm(const Vector& v)
{
	double max = 0, tmp = 0;
	for (size_t i = 0; i < v.size(); ++i) {
		tmp = fabs(v[i]);
		if (tmp > max)
			max = tmp;
	}
	return max;
}

double SquareNorm(const Vector& v)
{
	double sum = 0;
	for (size_t i = 0; i < v.size(); ++i)
		sum += v[i] * v[i];
	return sqrt(sum);
}

SymmetricTridiagonal::SymmetricTridiagonal(int _n) : n(_n) {}

SymmetricTridiagonal::SymmetricTridiagonal(int _n, double d, double a) : n(_n), diag(Vector(d, n - 1)), alpha(Vector(a, n - 2)) {}

SymmetricTridiagonal::SymmetricTridiagonal(int _n, const Vector& d, const Vector& a) : n(_n), diag(d), alpha(a) {}

const Vector SymmetricTridiagonal::operator*(Vector& x) const
{
	if (x.size() != n - 1)
		throw std::length_error("SymmetricTridiagonal: sizeof(x) != n-1");
	Vector y(n - 1);
	if (n <= 2)
		y[0] = x[0] * diag[0];
	else {
		y[0] = x[0] * diag[0] + x[1] * alpha[0];
		y[n - 2] = x[n - 3] * alpha[n - 3] + x[n - 2] * diag[n - 2];
		for (int i = 1; i <= n - 3; ++i)
			y[i] = x[i - 1] * alpha[i - 1] + x[i] * diag[i] + x[i + 1] * alpha[i];
	}
	return y;
}

const Vector& SymmetricTridiagonal::GetDiag() const
{
	return diag;
}

const Vector& SymmetricTridiagonal::GetAlpha() const
{
	return alpha;
}

PossionMatrix::PossionMatrix(int _n) : SymmetricTridiagonal(_n, 2.0 * _n * _n, -1.0 * _n * _n) {}



RHS::RHS(int _n, double left_boundary, double right_boundary, std::function<double(double x)> func) : n(_n), v(Vector(n - 1))
{
	for (int i = 0; i < n - 1; ++i)
		v[i] = func((1.0 + i) * 1.0 / n);
	v[0] += left_boundary * n * n;
	v[n - 2] += right_boundary * n * n;
}

const Vector& RHS::GetVector() const
{
	return v;
}

const int RHS::GetSize() const
{
	return n;
}

double& RHS::operator[](std::size_t i)
{
	return v[i];
}


WeightedJacobi::WeightedJacobi(const PossionMatrix& A, const Vector& _f, double _omega) : n(_f.size() + 1),
	T(SymmetricTridiagonal(n, 1 - omega / 2.0 / n / n * A.GetDiag(), -omega / 2.0 / n / n * A.GetAlpha())),
	c(omega / A.GetDiag() * _f), omega(_omega) {}

void WeightedJacobi::IterateOnce(Vector& v) const
{
	v = T * v + c;
}

std::ostream& operator<<(std::ostream& os, const Vector& v)
{
	os << fixed << setprecision(3);
	os << "(";
	for (size_t i = 0; i < v.size() - 1; ++i) {
		os << v[i] << ", ";
		if (i % 10 == 0 && i > 0)
			os << std::endl;
	}
	os << v[v.size() - 1] << ")";

	return os;
}
